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Abstract 



The subsystem formulation of linear response time-dependent DFT 
(TD-DFT) is analyzed in detail. Two equivalent derivations are pre- 
sented and naturally yield self consistent subsystem TD-DFT equa- 
tions. One derivation is equivalent to the subsystem TD-DFT for- 
malism of Neugebauer [J. Neugebauer, J. Chem. Phys. 126, 134116 
(2007)]. The other derivation yields Dyson type equations involving 
three types of subsystem response functions: coupled, uncoupled and 
Kohn-Sham. 

Analysis of the pole structure of the subsystem response functions 
shows that a single subsystem response function contains information 
about the electronic spectrum of the entire supersystem. Compari- 
son of the subsystem and supersystem response functions shows that, 
while the correlated response is subsystem additive, the Kohn-Sham 
response is not. This particularly puzzeling aspect of the theory is 
found to be an artifact stemming from the subjective nature of the 
density partitioning in subsystem DFT. This aspect is analyzed and 
compared to the less subjective partition DFT theory. 

This work sets the stage for a number of theoretical advances, 
such as the development of quasiparticle GW corrections to subsystem 
DFT, the development of a mapping between subsystem and supersys- 
tem Kohn-Sham orbitals, and the evaluation of non-local correlation 
energy with the adiabatic connection fluctuation dissipation theorem. 
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1 Introduction 



When modeling systems that contain a large number of electrons, even the 
Kohn-Sham Density Functional Theory (KS-DFT) approach [1] has its hm- 
its. In the past decade many approximations, such the Fast Multiple Method 
of Scuseria [2] or methods exploiting the locahty of chemical bonds [3] , as well 
as the so-called real-space methods [4] made it possible to massively reduce 
KS-DFT complexity for spatially extended molecules. However, the large 
pre-factor of such scaling laws left the calculation of most realistic, fuUy- 
solvated systems still prohibitive [5,6]. Despite the relatively favourable 
computational scahng of KS-DFT with the system size, most of KS-DFT 
simulations are bound to consider only small portions of the larger and more 
realistic system. 

Reducing the computational complexity of KS-DFT by partitioning the 
total electron density of a system into subsystem contributions has been an 
appealing idea since the early works of Gordon and Kim [7, 8] . However, 
the success of KS-DFT seemed to have rendered partitioning methods un- 
necessary. This is evident from the Quantum Chemistry literature of the 70s 
and 80s, where partitioning methods were frequent only to high-end wave 
function methods, and interactions between subsystems were treated with 
various types of perturbation theory [9] . Despite two successful applications 
of density partitioning techniques, first by Senatore and Subbaswamy [10], 
and then by Cortona [11], revival of density partitioning methods is due to 
a paper by Wesolowski and Warshel published in 1993 [12]. Presently, sub- 
system DFT is being developed by many research groups worldwide [13-21]. 
Successful apphcations of subsystem DFT are reported for applications re- 
lated to the ground state, such as analysis of electron densities [22] , and spin 
densities [23]; and for calculations of charge and excitation energy transfer 
parameters [15,24,25]; as well as for electronic spectra and molecular prop- 
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erties [20,26-30]. 

The time-dependent extension of susbsystem DFT has been pioneered 
by Casida and Wesolowski [31]. However, Neugebauer [32] is credited for 
deriving working equations for the solution of the subsystem time-dependent 
DFT (subsystem TD-DFT, hereafter) and for applying subsystem TD-DFT 
to determine excitation energies [20,28,32], charge/exciton couplings [15,24, 
25] , and molecular properties [33] . Similarly to subsystem DFT, subsystem 
TD-DFT is developed to take full advantage of the subsystem nature of 
the majority of real life systems. Solvated systems are a typical example 
of this. An early success story of subsystem TD-DFT is the calculation 
of solvatochromic shifts [34]. More recently, the electronic spectra of light 
harvesting complexes model systems containing more that 1000 atoms has 
been calculated with susbsystem TD-DFT [28]. 

Such a large body of work on subsystem DFT and TD-DFT shows their 
usefulness and importance. Therefore, rather than focussing on proving that 
susbsytem TD-DFT is a valuable alternative to regular TD-DFT of the su- 
persystem, this work aims at clarifying the relationship between the subsys- 
tem and regular TD-DFT of the supcrsystem. Specifically, two questions are 
answered: (1) How docs the KS and interacting response functions of the 
supcrsystem relate to the KS and interacting response functions of the sub- 
systems? And (2) How do practical subsystem TD-DFT schemes compare 
to the exact theory? 

Once a relationship between the subsystem formulation and the TD-DFT 
of the supersystem is clarified and exemplified by rigorous equations, a num- 
ber of important advances of subsystem DFT will become available. For 
example, it will be possible to apply subsystem TD-DFT to compute non- 
local correlation energy with the adiabatic connection fluctuation dissipation 
theorem [35,36], and to develop quasiparticle perturbativc corrections to sub- 
system DFT, similarly to what done by the GW approximation for regular 
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KS-DFT [37]. 

In this work, a simple and general formalism of subsystem TD-DFT is de- 
rived, including equations needed for practical calculations involving pertur- 
bative solutions. The new derivations will be amenable to a deeper analysis 
and understanding of the theory of subsystem TD-DFT. In addition, rela- 
tion with the TD-DFT of the supersystem will be established and analyzed. 
Dyson type equations for susbsystem TD-DFT are sought and are used to 
determine how approximate treatments compare with the exact formalism. 

This work is organized as follows. In the next section, a theoretical back- 
ground is given on KS-DFT, subsystem DFT, and hnear response TD-DFT. 
In Section 3, a rigorous derivation of linear response subsystem TD-DFT is 
carried out. In Section 4 an alternative derivation of subsystem TD-DFT is 
presented in terms of subsystem response functions. Section 5 is devoted to 
the comparison of subsystem TD-DFT with TD-DFT of the supersystem. In 
Section 6, conclusions and future directions are drawn. 

2 Theoretical Background 

2.1 Ground state DFT and subsystem DFT 

KS-DFT can be summarized by the following equation, known as the KS 
equation, in canonical form, 

0fc(r) = Sk<pkir), (1) 

where Vcs is the effective potential that the one-particle KS orbitals, (f)k, 
experience, and are the KS orbital energies. The spin labels have been 
omitted for sake of clarity, as throughout this work only the spin restricted 
case is considered without loss of generality of the derivations. The electron 
density is simply p(r) ^2J2T'' 
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The effective potential, Ves, is given by 



Ves{r) = Ve^t{r) + WeN(r) + i;coui(r) + v^^ir), (2) 



with Vext being an externally applied time-independent potential, the 
electron-nucleus attraction potential, "i^coui the Hartree potential, and fxc the 
exchange-correlation (XC) potential. Usually [1] , is part of f ext as they 
are both density independent. However, for sake of a consistent presentation 
of the following derivations, the two potentials are kept separated. 

Subsystem DFT is based on the idea that an electronic (molecular) sys- 
tem can be more easily approached if it is partitioned into many smaller 
subsystems. In mathematical terms, this is done by partitioning the electron 
density as follows [10, 11] 

Ns 

p{r) = ^p,(r), (3) 

I 

with Ns being the total number of subsystems. Self consistent solution of 
the following coupled KS-like equations (also called KS equations with con- 
strained electron density [38]) yield the single subsystem KS orbitals, i.e. 



<f>i{r)^ei<f>i{r), with / = 1, . . . , TV^ (4) 
with the effective subsystem potential given by 

^^eff(r) = 'yext(r) + vi^{r) + t;coui(r) + vi^{r) +vl^y,{r) . (5) 



same as regular KS— DFT 



In the above it is clear that if an external potential is applied to the system, 
Vext, every subsystem will experience the same potential. In the so-called 
Frozen Density Embedding (FDE) formulation of subsystem DFT [12,38], 
the unknown potential above, Vemh, is called embedding potential and is 
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given by 



^emblr 



5p{r) 5pi{r) 5p{r) 5pi{r) 



(6) 



Throughout this work, "subsystem DFT" is used as a synonym of FDE. 
The density of the supersystem is thus given by p(r) = 2^^^^ YIT'^' |0f (^) f- 

2.2 Linear response TD-DFT 

The time-dependent KS equation, 



(r,i)-z- 



(r,t) 



dt 



(7) 



relates the time dependent KS orbitals, ^^(r, i), and the correlated density, 
p{r,t) — 2 J^^'^^ |0fe(r, with the external perturbation Uext(i", When 
starting from the ground state density, the time-dependent KS potential is 
defined according to Eq. (2) as 



Ves{r, t) = Vext(r, t) + VeN(r) + i;coui(r, t) + Vxc(r, t), 



(8) 



The external time-dependent potential in Eq. (8) is often referred to as "ap- 
plied" potential [31,39] and it constitutes the only time-dependent perturba- 
tion causing the density p to become a time-dependent function. With the 
exception of Wcn, which is considered static (the nuclei are assumed to be 
still in the time the perturbation is applied), the other potential terms part 
of the effective KS potential are dependent on time, but only as a result of 
the perturbation t;ext(i', 

Linear-response TD-DFT is based on the assumption that the density 
response to the external weak perturbation is given by the following linear- 
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response integral equations [40] 




(10) 



(9) 



where 



5Veff(r', t) = 5vext(r', t) + 5Vind{r\ t), 



(11) 



The induced potential, Svi^d, is expressed in linear response as well, namely 



The quantity ^g^l^}^ is called the XC kernel, /xc(r', t, r", t'). The functions 
x(r,t,r',t') and x^lr, r', t') are the correlated and the simplified KS re- 
sponse functions (or simply "correlated response" and "KS response", re- 
spectively). Eq. (9) constitutes the definition of linear- response TD-DFT, 
and Eq. (10) derives from Eq. (9) from the Runge-Gross theorem (Theorem 
4 of Ref. [41]). In the rest of this work, retardation effects are neglected and 
the adiabatic approximation is assumed, i.e. /xc(i", t, r', t') — /xc(r, r')5{t — t'), 



X^{r,t,r',t') — x°(r, r', t)5(i — t'). In addition, 5vext{i^,t) — ^fext(r) in the 



linear response integral equations, and thus the XC kernel is assumed to be 
exactly the functional derivative of the XC potential [42]. 

For practical calculations, Eq. (10) is the most important. This is because, 
in the adiabatic approximation it involves quantities that can be easily ex- 
tracted from the ground state KS system, such as the KS response function 
(given here in Fourier transform) 



Jt'=to J Ll"^^ ~ 




5p{v",t')dv"dt'. 



(12) 



occ virt 



X°(r,r',a;) = J] J] V^'^»(r)0a(r)0,(rO0a(r'), 



(13) 



i a 
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where 0i and (j)a are occupied and virtual KS orbitals. In a simpler notation, 
omitting the integral signs and the variable dependence, practical calcula- 
tion of the density response are carried out by self-consistently solving the 
following [40] 

\xr'-f[ 

with 

/(r,r') = ^^ + /xc(r,r'). (15) 
As Eq. (14) must hold for any Svext, comparison with Eq.(9) yields 

(xr' = (16) 

also known as the Dyson equation for the response function [40,42,43]. 

3 Subsystem TD-DFT 

This section is devoted to the derivation of linear response subsystem TD- 
DFT. Even though this theory has been first derived by Neugebauer [32], 
here it is presented in a different mathematical formalism. The derivations 
and analyses presented in this section are important as they pave the road 
to the formalism presented the subsequent sections. 

3.1 Mathematical derivation 

Following the usual decomposition of the density change in subststem TD- 
DFT [20,25,32], the total electron density change of the system, 5p, due to 
an external perturbation is given exactly by 

Ns 

6p{r,t) = J2^Pi{r,t), (17) 
I 



(14) 
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where 5pi is the density change of the single subsystem I. Similarly to 
Eq. (11), consider the effective time-dependent perturbation on subsystem /, 
^^^^(r,^), being a functional of all the subsystem densities, and defined as 
follows 

5vi^{v,t)^5v^Jj^+ 5vUr^ . (18) 

perturbation . , , , , 
induced potential 

on subsystem I 

Similarly to Eq. (5) , in the above it is assumed that the external potential is 
applied to the entire system and therefore it is the same external potential, 
5vcxt, that interacts with all the subsystems. The induced potential, 5^^4d' 
be defined in terms of the following functional derivatives of the subsystem 
KS potential given in Eq. (5) [20] 

^-X + ^-*l(r,r') + ^^-^4f^L^„. (20) 

dpj{v') |r — r'l dp{r)dp{r) dpi{r)opj{r') 

Following the adiabatic approximation, in the above equation the definition 
of XC kernel, /xc from Eq. (12) has been used and the time-dependence of 
the various kernels (kinetic and XC) has been neglected. Similarly to Eq. (9) 

and Eq. (10), with the aid of the Runge-Gross theorem, the time-dependent 
subsystem density can be obtained self consistently as 

6pi{r,t) = J x'i{r,r',t)6vUr')dr' (21) 
= J x?(r,r',t)H'ff(r')dr', (22) 

where x? is the KS response of the subsystem to the external perturbation, 

Xj is the correlated "coupled" subsystem response function, and Svl^ is given 
by Eqs. (18-19). 
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The above equations hold a great deal of information, e.g. the subsystem 
time-dependent density can be obtained from the simplified subsystem KS 
response function and the effective time-dependent potential. Eqs. (18-19) 
can be used in Eq. (22), yielding 

Spi{r,t)= J x%r,r',t)Sv,^t{r')dr'+ 



and by grouping the terms involving Spi on the Ihs and expressing 6pi{r,t) 
in terms of integrals of suitable Dirac deltas 



5(r-r')5(r'-r")-X?(r,r',t) 



5pi(r",t)dr'dr" ^ 



J x?(r,r',i)5^ext(r')dr' + J ;^?(r, r', t) ^^^p,(r", i)dr'dr 



(24) 



Eq. (24) can be rearranged by acting on the left by (x5) ^ (r'''^ r, t) and inte- 
grating over dr, using the relation J {x^)~^ (r'", r, t)x%r, r', t)dr = 5{r"' — r'), 



J [(X?)"' (r'",r',t)5(r' - r") - ^(r'" - r')^" 
= 5^ext(r'") + 1 6{r"' - r') J] ^-^^Spj{r", t)dr'dr". (25) 

After integration over r' and substitution of r'" — > r and r" — > r' the following 
is obtained 



5pi{r",t)dr'dr" 



J SMr',t)dr'^Sv..,{r)+jY^'-^Spjir',t)dr' 



(26) 
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Substitution of 



allows for the definition of the "uncoupled" subsystem response function as 



(x^)-^(r,r',t) = {x':)-'{ry,t)-Kn{ry). 



(28) 



Reahzing that Eq. (26) holds for every subsystem, the following Ng x Ns 
matrix vector equation can be formally constructed 



with 



M 



and 



MSp = iSv^^t, 

( \ 



(29) 



(30) 



(31) 



\ ^Pns I 

where K is the matrix composed of the Ku kernels. If the matrix in Eq. (30) 
is invertible, then the poles of occur at the true excitation energies of 
each subsystem, and hence of the total supersystem. Eq. (29) can be con- 
sidered the subsystem DFT equivalent of Eq. (14). The matrix formulation 
above yields the coupled subsystem response function defined in Eq. (22) as 



(32) 



thus a formal relationship is 



(33) 
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and 

5p = Tr [M-i] 5v,,f (34) 

Equations (32-34) are well suited to be used in practical calculations. This 
is because, in practice, all the operators (response functions and kernels) are 
expressed in a matrix form. The next section is devoted to the exploration 
of practical schemes for the solutions of Eqs. (29-34). 

3.2 Practical calculations 

In the occupied . . .) and virtual (a, b, . . .) KS orbital basis set, the su- 
permatrix K occurring in Eq. (30) is defined as 

^i^aMJ»h = J </'..(r)0a,(r) (j^) <PjA^')cP,A^')dr dr'. (35) 
With that, the orbital rotation Hessians can also be defined in a finite basis 
(A - B)(^ia)j(jb)j = {ei - el) 6ij6abSij, (36) 

(A + B)^ia)j(jb)j = {ei - €•) SijSabSlJ + 2K^ia)j(jb)j- (37) 

If X and Y are vectors associated with the usual excitation and de-excitation 
components of the transition density matrix [44,45] , solutions of the following 
anti-hermitian eigenvalue problem [28] corresponds to setting (^Vext — and is 
equivalent to inverting the matrix in Eq. (30) , for only two subsystems using 
Eq. (29) yields 
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where, in the finite basis, dpj — \ , A/ = , and A/ are 

Y, ; 1^ B, A, 

the matrices with zero off-diagonal element and 1 in the diagonal elements 
corresponding to X/ and —1 for the ones corresponding to Y/. A perturbative 
approximation to Eq. (33) is [25,32] 



(^[Ai - 00 Ai] - K/j [A^ - ouAj]-' Kjj^ Spj = 0. 



(39) 



The matrix representation of the subsystem response function in the un- 
coupled formalism is defined as 

X'i = [Aj-u;Aj]-\ (40) 

Applying Eq. (39) itcratively yields the subsystem coupled response {xj) as 
defined in Eq. (33). If, instead, it is applied only once, it yields a perturbative 
subsystem coupled response function 

Ns 

fepert)^' = Uir' -Y.^IjX'-jKji. (41) 

The use of Eq. (40) and iteratively of Eq. (39) [or equivalently Eq. (41)] for 
the solution of uncoupled and coupled subsystem TD-DFT, respectively, has 
shown great success for the calculation of excitation energies, as well as for 
such molecular properties notoriously affected by the environmental response, 
such as circular dichroism and excitonic couplings [20,25,27,28,32,33,46]. 

Several remarks should be made about the above formalism. In the de- 
generate case, for example when dimers are considered, the above equation 
runs into numerical inaccuracies. That is because of the occurrence of poles 
in the uncoupled response at exactly the same frequencies as the zeroes of the 
inverse uncoupled response. In those cases, Neugebauer has circumvented the 
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problem [32] by abandoning the effective perturbative treatment above and 
by directly inverting the part of the M matrix corresponding to the subspace 
spanned by the two monomers. In a second step, Eq. (41) can be applied to 
account for the coupling with all the remainder subsystems. In the degen- 
erate case, the coupled subsystem response function is not given by Eq. (39) 
but instead takes a more complicated form, given implicitly by Eq. (34) . It 
is, therefore, concluded that Eq.(41) can be considered a perturbative solu- 
tion [47] and thus it does not represent the true coupled subsystem response 
function. 

In the following section, an alternative derivation of subsystem TD-DFT 
is carried out. The new derivation is such that: (1) the coupled response is 
defined explicitly even in the degenerate case, and (2) it is consistent with 
and recovers the original formulation. 

4 Alternative Derivation of Susbsystem TD- 
DFT 

4.1 Mathematical derivation 

Eq. (26) can be rearranged as follows, 

J (x^)-' (r, r', t)xKr', r", t)5v,^,{r")dr"dv' = 

5(r - r')5(r' - r") + ^ Kjj{r, r')x}(r', r", t) 5^;ext(r")dr"dr'. (42) 
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The above equation must hold for any (5vext(i""), and specifically for 5vext(r") = 
5{r" — f). Integration over dr" yields 




(43) 

Thus, after applying (r", r, t) and integration over dr, the following Dyson- 
type equation is obtained, in simplified notation, 

Ns 



The main difference between the above Dyson equation and the one in 
Eq. (41) is that the former should be considered a more general Dyson equa- 
tion, as it is derived in the general case without assuming a finite basis expan- 
sion of the involved operators. In addition, it does not employ a perturbative 
approach. 

Dyson equations for the response functions involving only the kernels 
and the KS response functions are derived starting from Eq. (26) and read as 
follows 

XI-XU x'Knxl (45) 

Ns 

X/ = X? + x'J2^ijXj- (46) 
J 

Similarly to regular TD-DFT, this formulation shows that the uncoupled 
response in Eq. (45) is similar to the one of the isolated subsystem, with a 
small correction in the kernel due to the second functional derivative of the 
kinetic energy. In addition, the coupling between subsystem responses in 
Eq. (46) is mediated by the off-diagonal elements of the kernel matrix K. 
In sharp contrast to the original formalism, here it is evident that if the 
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poles of the response function of subsystem / are well separated from the ones 
of the other subsystem response functions, then the poles of each subsystem 
response contain the ones of all other subsystems. This is a particularly 
interesting result, as it shows that formally the correlated response function 
of a single subsystem contains information about the electronic spectrum 
of the entire supersystem. Obviously, in the limit of infinitely separated 
subsystems, Kjj {r,r') will be identically zero when r and r' span regions of 
space occupied by different subsystems. Thus, the above observation needs 
to be taken with a grain of salt as it is valid only if the subsystems are 
relatively close to each other. 

Another interesting outcome of this formalism is that when two subsys- 
tems have poles at the same frequencies in the isolated case (or in the uncou- 
pled case) , then this degeneracy must disappear in the coupled case otherwise 
the response function would feature an unphysical "double pole". This im- 
plies that the above formalism is coherent with the existence of Davydoff 
splittings in dimeric systems [32,48]. 

4.2 Practical calculations 

Even though Eq. (44) is aesthetically pleasing, it is not suitable for practical 
calculations. The perturbative solution in Eq. (29) is recovered by rewriting 
the above equation in matrix form. In the first step, let us rewrite Eq. (44) 
as 



In the second step, the above equation is written in a more convenient matrix 
form 




(47) 



1 = Mx^ 



(48) 
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where 



-1 



(49) 



is the same matrix defined in Eq. (30) . Similarly as before, from Eq. (48) , the 
poles of are also the poles of the coupled response function. With this, 
the above matrix inversion can be approximated by a perturbative approach 
similarly to what is done in Eq. (41). 



5 Comparison to TD-DFT of the supersys- 



Comparison of subsystem DFT with regular KS-DFT is straightforward. In 
subsystem DFT one has to solve coupled KS-like equations, where the cou- 
pling term is conveniently expressed as a potential term, Vemh, added to the 
KS effective potential of the isolated subsystem. This means that the two 
formalisms involve similar algorithms for practical calculations. 

As it will be clear from the following derivations, this is not the case 
for the time-dependent extensions. In the following, Dyson-type equations 
will make it possible to directly compare subsystem DFT with the TD-DFT 
of the supersystem. The correlated response of the super system has a sim- 
ple relationship to the subsystem correlated responses, taking the functional 
derivative with respect to 6vext of both sides of Eq. (17) one obtains 



Conversely, due to the non- uniqueness of the density partitioning in Eq. (3), 
the "simplified" KS response of the supersystem has no simple relationship 



tern 




(50) 



I 
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with the subsystem KS responses. 



5.1 Subsystem versus full KS response function 

In order to find a relationship between the subsystem and the supersystem 
KS response functions, let us manipulate Eq. (21) by inverting the subsystem 
response functions, one by one, in the following manner 

(x?)"'5p7 = H'ff, (51) 

where we have omitted the integration symbols for sake of a lighter notation. 
An important difference between the induced potential used in subsystem 
TD-DFT defined in Eq. (20) and the one used in the TD-DFT of the super- 
system defined in Eq. (12) resides in the kinetic energy kernels. The two can 
be related, 

Svif,{r, t) = 5v,s{r, t) + ^^(r, t), (52) 

where, using Eq. (19) and Eq. (20), we define the kinetic energy part of the 
subsystem kernel as 

J J V^P(r)^p(r') 5pi{y)5pj{y')J 
= j Mr,r')5p{r',t)- J fi{v,v')5pi{v' ,t), (53) 

where Eq.(17) has been used for the first term of the rhs. Using Eq.(53) in 
Eq. (51), we obtain 



(X?)"' + Spj = (1 + /tX°) Sv^S: (54) 



where the number 1 above is intended to be the identity in functional space, 
i.e. a Dirac delta in the position representation. Inverting the operator on the 
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Ihs of the above equation and summing over all the subsystems, we obtain 
the following 



11^ ^ 



(55) 



At this point there are several algebraically non-equivalent ways to proceed. 
Two routes are considered here: the first one leading to an exact expression, 
and the second one leading to expressions suited for approximations. Section 
5.2 is devoted to an analysis of the following derivations. 



5.1.1 Exact expression 

Extracting x° from the braces of Eq. (55), and realizing that x^^Vefr = 5p 



r r n -1 r n ^ 



By defining (xj) ' = [(x?) ' + /^l , Eq.(56) leads to 



eff- 



(56) 



X 



Xi) 



(57) 



An important remark is that the KS supersystem considered here is the true 
KS supersystem. Other formulations of subsystem TD-DFT [31], instead, 
considered a supersystem defined similarly to KS-DFT, but solved with a 
modified potential containing the KS potential plus the total kinetic energy 
potential, 
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5.1.2 Approximate expressions 

A first approximation can be reached directly from Eq. (57) in the Hmit of 
vanishing kernels, namely 

Ns 

(X°)-' = E(X?)"- (58) 

However, a different approximation can be make by first taking the func- 
tional derivative with respect to Sves on both sides of Eq. (55), namely 

Ns 

/ = Ex?[1 + /^X?]"'(1 + /t/), (59) 
which can be arranged to 

= [I - frS]-' S, (60) 

with S = X^/'^X? [l + /tX?] ^- '^^^ above inverse operations expression can 
be approximated with linear expansions, in the limit of small /fx? small 
/tX?, to 

Ns Ns Ns 

/ ^ E^? - + E x'Jtx'j, (61) 

II IJ 

featuring an interesting resemblance to the Dyson equation for the response 

function. 

It is instructive to notice how the supersystem KS response function cal- 
culated with the linear expansion of the inverse operations in Eq. (61) and 
neglecting the kinetic and exchange-correlation kernels, which we term here 
Random Phase Approximation (RPA), is exactly the sum of the subsystem 
response functions, i.e. 

Ns 

^0,RPA ^ ^^0,RPA (g2) 
/ 
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Assessing the performance of the RPA approximation compared to the exact 
expression in Eq. (57) and the approximation in Eq. (58) hes beyond the scope 
of this theoretical work and will be the focus of a future work. 

5.2 Physical meaning of the subsystem KS responses 
and comparison to PDFT 

The derivations in the preceding section stand out as being too complicated 
for just the KS response, paradoxically in this context known as the "simpli- 
fied" response. What is the significance of such a comphcated relationship 
between the supersystem and the subsystem KS response functions? What 
is puzzehng in Eqs. (57-61) is that the KS response of the supersystem con- 
tains terms coupling the subsystem KS responses. This is not a very good 
property of this theory, as subsystem additivity is sought in the density, in 
the correlated response in Eq. (50) and it is expected to appear in the KS 
density response as well. 

This apparent artifact is due to the non-uniqueness and subjectivity of 
the density partitioning employed in Eq. (3). An indication of this artificial 
behavior of the subsystem KS responses can be easily shown by considering 
a more refined version of subsystem DFT known as partition DFT (PDFT). 
In PDFT theory [49-52] the effective subsystem time-dependent potential is 



where 5v^ is the change in the partition potential (a quantity shared by 
all subsystems and thus unique). The above equation can be rearranged 
similarly to the step carried out between Eq. (54) and Eq. (55) , to yield 



(63) 



(64) 
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with fp = The above equation is rearranged to 



1 



(65) 



which can be approximated assuming small fpXi by 



Ns Ns 



(66) 



/ I 



The last two equations feature no cross terms coupling the subsystem KS 
responses. Thus, PDFT provides a more intuitive time-dependent behavior 
of the subsystems and is completely free of artifacts due to the non-unique 
partitioning appearing in regular subsystem DFT. 

Non-orthogonahty also plays a role. For example, the simple relationship 
recovered in the RPA approximation in Eq. (62) is expected to hold more 
strongly in the limit of small electron density overlap between subsystems, 
i.e. quasi orthogonality, as the non-additive functionals in that case are close 
to being identically zero. 

Another important remark is that Eqs. (57-61) establish a map between 
KS and FDE orbitals without having to invert the KS equation as it is done 
in recent works [53,54]. In principle, Eq. (61) can be used to approximate the 
KS orbitals by correcting to first order the subsystem orbitals. Developing 
the working equations for such a mapping is not in the scopes of this work 
and will be the subject of future studies. 

6 Conclusions 

In this work, the theory of hnear-response subsystem TD-DFT is analyzed 
in detail. Two formulations are considered, one is equivalent to the theory 
presented by Neugebauer [25] , and the other is presented in terms of Dyson 
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equations involving subsystem response functions. Three types of subsystem 
response functions are considered: coupled, uncoupled and KS response. The 
coupled and uncoupled are exact and approximated correlated subsystem 
responses, respectively. 

It is found that, for non-infinitely separated subsystems, the pole struc- 
ture of the correlated (coupled) subsystem response functions contains the 
excitations of the entire supersystem. This shows that if an external potential 
is in resonance with an electronic transition of one subsystem, the electronic 
response of another subsystem will also be strongly affected. This behaviour 
generally does not fit a picture of "locahzed excitations" [16] but instead is 
consistent with the idea that the response of a collection of subsystems is col- 
lective. Localization of the excitations may take place whenever the kernel 
coupling the subsystem's excitations Ku is small. 

Comparison of the subsystem responses with the response functions of the 
supersystem shows that, while the correlated response is subsystem additive, 
the KS response is not. This partucularly puzzeling aspect of the theory is 
found to be an artifact stemming from the subjective nature of the density 
partitioning, as it is not found in time-dependent PDFT theory. 

This work sets the stage for a number of important theoretical advances, 
such as the development of quasiparticle GW corrections to subsystem DFT, 
and the setting up of a mapping between subsystem and supersystem KS 
orbitals through the expressions relating the KS response functions. The 
evaluation of non-local correlation energy, and thus the account of dispersion 
interactions between subsystems, is also a possible outcome of this theory 
and is being investigated as a future development. 
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